AEI-2007-140 



Type II critical phenomena of neutron star collapse 
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We investigate spherically-symmetric, general relativistic systems of collapsing perfect fluid distri- 
butions. We consider neutron star models that are driven to collapse by the addition of an initially 
"in-going" velocity profile to the nominally static star solution. The neutron star models we use are 
Tolman-Oppenheimer-Volkoff solutions with an initially isentropic, gamma-law equation of state. 
The initial values of 1) the amplitude of the velocity profile, and 2) the central density of the star, 
span a parameter space, and we focus only on that region that gives rise to Type II critical behavior, 
wherein black holes of arbitrarily small mass can be formed. In contrast to previously published 
work, we find that — for a specific value of the adiabatic index (F = 2) — the observed Type II critical 
solution has approximately the same scaling exponent as that calculated for an ultrarelativistic fluid 
of the same index. Further, we find that the critical solution computed using the ideal-gas equations 
of state asymptotes to the ultrarelativistic critical solution. 
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I. INTRODUCTION 



Critical phenomena in general relativity involves the 
study of solutions — called critical solutions — that lie at 
the boundary between black hole-forming and black hole- 
lacking spacetimes. (See P, H, Q for reviews.) Published 
work in general relativistic critical phenomena began over 
a decade ago with a detailed numerical examination of 
the collapse dynamics of a massless scalar field, minimally 
coupled to the general relativistic gravitational field [J]. 
This first study in critical phenomena touched upon the 
three fundamental aspects of black-hole-threshold criti- 
cal behavior: 1) universality and 2) scale invariance of 
the critical solution with 3) power-law behavior in its 
vicinity. All three of these features have now been seen 
in a multitude of matter models, such as perfect fluids 
@, S 0, an SU(2) Yang-Mills model @, [|, and colli- 
sionless matter [13, EI to name a few. It was eventually 
found that there arc two related yet distinct types of crit- 
ical phenomena: Type I and Type II, so named because 
of the similarities between critical phenomena in general 
relativity and those of statistical mechanics. 

Type II behavior was the first to be discovered 
and entails critical solutions that are either continuously 
self-similar (CSS) or discretely self-similar (DSS). Super- 
critical solutions — those that form black holes — give rise 
to black holes with masses that scale as a power-law, 
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implying that arbitrarily small black holes can be formed. 
Here, p parameterizes a 1-parameter family of initial data 
with which one can tune toward the critical solution, lo- 
cated at p = p* , and 7 is the scaling exponent of the 
critical behavior. Since Mbh(p) — > as p — > p* , this 
type of critical behavior was named "Type II" since it 
parallels Type II (continuous) phase transitions of statis- 
tical mechanics. 

As in the statistical mechanical case, there is also a 
Type I behavior in gravitational collapse, where the black 
hole mass "turns on" at a finite value. As might be ex- 
pected, Type I critical solutions are quite different from 
their Type II counterparts, tending to be meta-stable 
star-like solutions that are static or periodic. In this pa- 
per, attention is restricted to Type II behavior; results 
from our study of Type I transitions in our model are 
reported in a separate paper [ll|. 

The accepted picture describing the scaling behavior 
seen in Type II critical collapse was suggested by Evans 
and Coleman [|| , who computed the critical solution for a 
radiation fluid (fluid pressure, P, and density, p, related 
by P — p/3) in two distinct ways. First, using a code 
that solved the full set of partial differential equations 
(PDEs) for the fluid and gravitational field, and by tun- 
ing an initial data parameter as sketched above, Evans 
and Coleman were able to compute a strong field solu- 
tion that sat at the threshold of black hole formation, as 
well as establish a scaling law of the form ((T|) . Further- 
more, the results of this numerical experiment provided 
compelling evidence that the threshold solution was con- 
tinuously self-similar. Second, by adopting the assump- 
tion of continuous self-similarity as an ansatz, Evans and 
Coleman reduced the set of PDEs governing their model 
to a set of ODEs, from which a precisely CSS solution 
was calculated. The solutions computed using these two 
completely different techniques were found to agree ex- 
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tremely well. Crucially, it was argued that the observed 
scaling behavior of the mass above threshold could be 
explained using linear perturbation theory about a back- 
ground given by the CSS critical solution. Such an anal- 
ysis was carried out by Koike et al. [13| (for the radiation 
fluid), who showed that the scaling exponent, 7, was the 
inverse of the Lyapunov exponent of the critical solution's 
single, unstable eigenmode. 

Subsequent work showed that 7 was not a truly uni- 
versal constant, but that its value could depend on the 
specifics of the matter model used. The first evidence 
for this non-universality in scaling behavior was given in 
concurrent works by Maison [TJ and Hara et al. [HI]. 
Using similar methods to those of [H, Ojj], they found 
that 7 for an "ultrarelativistic" fluid with equation of 
state (EOS) 

p = (r-i) P (2) 

is dependent on the adiabatic index, T. 

Most of these investigations, however, have involved 
ultrarelativistic fluids that are explicitly scale-free. The 
reason for the predominance of this type of fluid is due to 
the fact that Cahill and Taub (lTJ showed that only those 
perfect fluids which have state equations of the form of 
Eq. @ — i.e. the so-called ultrarelativistic EOS — can give 
rise to spacetimes that admit a homothetic symmetry 
(i.e. which are self-similar). Hence, it is not completely 
unreasonable to expect that Type II, CSS critical solu- 
tions would only appear in such fluids, or at least in fluids 
that admit an ultrarelativistic limit. To study this con- 
jecture, Neilsen and one of the current authors [7j consid- 
ered the evolution of a typical perfect fluid with equation 
of state, 

P=(T-l)p e (3) 

that introduces a length scale into the field equations. 
Here, P is the pressure, p Q is the rest-mass energy den- 
sity, and e is the specific internal energy density. It was 
argued in Q that Type II critical collapse scenarios are 
typically kinetic energy dominated, entailing increasingly 
large central pressures that maintain the tenuous bal- 
ance between the matter dispersing from the origin and 
collapsing to a black hole. Therefore, if was reasoned 
that should one be able to give the fluid sufficient kinetic 
energy, then it would naturally enter an ultrarelativis- 
tic phase. Specifically, if the fluid undergoes a collapse 
such that e — * 00 dynamically, then p will effectively 
become negligible in the equations of motion (EOM) and 
the system will be able to follow a scale-free — hence self- 
similar — evolution. To see if this hypothesis was correct, 
compact distributions of perfect fluid, with P = 0.4p o e 
(r = 1.4) were collapsed, and the calculations were tuned 
to a threshold solution. The critical solution thus ob- 
tained by solving the full set of PDEs closely matched 
the precisely self-similar solution, which was calculated 
by assuming that a model governed by the ultrarelativis- 
tic EOM had an exact homothetic symmetry. Further, 



it was found that the scaling exponent, 7, defined by 
Eq. (TIJ matched that of the ultrarelativistic critical solu- 
tion with r = 1.4. Since the ultrarelativistic fluid exhib- 
ited Type II phenomena for all considered values of the 
adiabatic index in the range 1.05 < T < 2, the results of 
suggested that the Type II ideal-gas critical solution 
for any T in that range should be the same as that for 
an ultrarelativistic fluid with the same T. 

This hypothesis is not without precedence, since sev- 
eral models have been found to exhibit DSS or CSS col- 
lapse, even when explicit length scales are present. For 
instance, one of us found Type II behavior for the case of 
a collapsing massive scalar field [l|[ — that is a scalar field 
with potential V(cj)) — ^m 2 (f) 2 — even though the model 
has an explicit length scale set by 1/m. The heuristic 
argument presented in [lH ] is that the potential term is 
naturally bounded since <f> itself is bounded in the critical 
regime, but that the kinetic term — □(/) — diverges in the 
critical limit. Hence, the kinetic term overwhelms the po- 
tential term and essentially makes the critical evolution 
scale- free. 

The single study exhibiting Type II behavior in perfect 
fluid collapse with an ideal gas EOS [7j remained unveri- 
fied until work by Novak [19( . To determine the possible 
range in masses of nascent black holes formed from stellar 
collapse, Novak performed a parameter space survey us- 
ing Tolman-Oppenheimer-Volkoff (TOV) solutions with 
r = 2, varying the overall amplitude, £/ amp , of an other- 
wise fixed initial coordinate velocity profile for the fluid 
in order to generate critical solutions. The Type II be- 
havior observed was quantified by fitting to the typical 
black hole mass scaling relation |T]) , where p was identi- 
fied with C/ amp . Significantly, Novak was able to observe 
scaling behavior even with a realistic equation of state 
formulated by Pons et al. [20| . This was somewhat sur- 
prising, since there were some expectations that Type II 
phenomena would not be observed when using realistic 
equations of state [3j. However, provided that the equa- 
tion of state admits an ultrarelativistic limit, as is ap- 
parently the case for Novak's calculations, the heuristic 
argument sketched above suggests that one should expect 
to see Type II transitions. 

Although Novak observed Type II behavior, he did not 
find the same scaling exponent as had been observed for 
the r = 2 ultrarelativistic fluid in the study described 
in [7j. In addition, he claimed that 7 was a function of 1) 
the central rest-mass density, p c , which (as described in 
the next section) parameterizes the initial star solution, 
and 2) the EOS used. He observed that the fit to Eq. ([T]) 
worsened as p c increased to that of the maximum mass 
solution, and that it eventually broke down completely. 
Specifically, he found for the ideal-gas EOS © 

7 * 0.52, (4) 
and when using the realistic EOS 

7 ~ 0.71. (5) 

These values are significantly different from the values 
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most recently calculated with the r = 2 ultrarelativistic 
fluid [5] using a variety of methods: 



7 



0.95 ±0.1. 



(6) 



Here we have taken the average of the three independent 
values calculated in [5[ and the quoted uncertainty is the 
standard deviation of those values, and does not include 
the systematic errors inherent in the distinct calculations. 

However, Novak clearly states in [l9| that his code was 
not designed to simulate the formation of very small black 
holes, and apparently was only able to tune to a preci- 
sion of \p — p*\/p — 10~ 3 . In this paper, we reexamine 
the Type II behavior in this particular system in order 
to check the claims of [19(, and obtain what we claim is 
an improved measurement of the scaling exponent. Dif- 
ferent families of initial data are used to demonstrate the 
universality of our computed critical solution. Also, we 
compare the critical configuration calculated from a near- 
threshold neutron star collapse to the critical solution ob- 
tained using an explicitly ultrarelativistic fluid. To more 
accurately study CSS behavior as the black hole thresh- 
old is approached, we employ mesh refinement techniques 
and non-uniform discretization. Further, we implement 
methods that improve the accuracy of the transforma- 
tion from so called conservative variables to primitive 
variables that is required in our numerical analysis. 

The remainder of this paper is structured as follows. 
Section |TT] introduces the model and equations used to 
describe our collapsing neutron stars. In Sec. IIII1 we 
describe the numerical methods used to solve the coupled 
fluid and gravitational PDEs, and include a discussion 
of an instability we observe near the critical threshold. 
In Sec. IIV1 we analyze the observed Type II behavior 
and compare it to previously published work. We then 
conclude in Sec. IVl with some closing remarks and notes 
on anticipated future work. 

Geometrized units, G = c = 1, are used throughout, 
and our tensor notation follows 12111. 



II. THEORETICAL MODEL 

As in many previous critical phenomena studies in 
spherical symmetry 0, H, S B 0; we employ the so- 
called polar-areal metric 



ds 2 
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where a is often referred to as the lapse function. We use 
the perfect fluid approximation for the matter model of 
our neutron stars, so the stress-energy tensor takes the 
form 



Tab = {p + P) u a ii b + Pg ab 



(8) 



Here, u a (r,t) is the 4- velocity of a given fluid element, 
P(r, t) is the isotropic pressure, p(r, t) = p a (1 + e) is the 
energy density, p (r, f ) is the rest-mass energy density, 
and e(r, t) is the specific internal energy. 



Modern computational methods that cast the hyper- 
bolic fluid equations of motion in conservation form, and 
that use information concerning the characteristics of the 
equations, have been used very successfully in the mod- 
elingof highly-relativistic flows with strong gravity (see 
[221. |23|. [241 . I25I , [26j for a small, but representative, se- 
lection of papers on this topic). Here we adopt the for- 
mulation used by Romero et al. [26|], with a change of 
variables similar to that performed in [2|| . The formula- 
tion described in (26| has been used extensively for treat- 
ing relativistic flows in the presence of strong gravitation 
fields, and appears to work quite well in such instances. 

The EOM for the fluid are derived from the local con- 
servation equations for energy and baryon number which 
are, respectively, 



V Q T a b = 0, 



V a (p o u a )=0. 



(9) 



(10) 



Rather than working directly with the components of the 
fluid 4-velocity, it is more useful to employ the radial 
component of the Eulerian — or physical — velocity of the 
fluid as measured by an Eulerian observer: 



v(r, t) 



(11) 



where = [u*, u r , 0, 0]. The associated "Lorentz gamma 
function" is defined by 



W(r,t) = au l . 

and satisfies the usual relation 

1 



W 2 = 



1-v 2 ' 



(12) 



(13) 



since the 4-velocity is time-like and unit-normalized, 
i.e. u^u^ = —1. Adopting a notation where bold face 
symbols denote state vectors, the fluid EOM can be cast 
in conservative form as 



1 



d t q+-d r (r 2 Xf) =ip, 



(14) 



where q = q(r, t) is a state vector of conserved variables, 
and f = f(q) and ip = i/>(q) are, respectively, the flux 
and source state vectors. Our choice for the conserved 
variables is the one used almost exclusively in the field 
U 0, HI US 13, namely q= [D(r,t),S(r,t),r(r,t)] with 



D = ap W, 

S = p hW 2 v, 

t = E-D, 

E{r,t) = Po hW 2 -P, 



(15) 



and where h(r, t) = l + e + Pj p is the specific enthalpy of 
the fluid. D, S, E, and r can be thought of as the rest- 
mass density, momentum density, total energy density, 
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and internal energy density, respectively, as measured in 
a Eulcrian-frame defined by the ADM slicing. We found 
that for extremely relativistic flows near the threshold of 
black hole formation, this formulation was not very sta- 
ble. We therefore use a different formulation motivated 
by [H| j where it was found that evolving r ± S allowed 
for a more precise calculation since r ~ S in the ultra- 
relativistic regime. We thus define new variables 

n(r,i) = t + S, $(r,*) = t-S (16) 

and the state vectors become 



</> = 



D 

n 
$ 


E 

-E 



Dv 

v(U + P) + P 
w($ + P) -P 

P 

v 

Po 



(17) 



The elements of the vector w are the set of primitive 
variables used. The source function E takes the form 



e 



2PX 



where 



e 



aa 



(Sv -E)( 



8nrP 



P- 



(18) 



(19) 



With this source, the governing equations with which 
we solve for our metric functions are the Hamiltonian 
constraint of the ADM [27J formulation 



a 
a 



1 

4nrE 

2r 



1 

27 



and the polar-areal slicing condition 



a 
a 



47rr (Sv + P) + -j 



(20) 



(21) 



In order to monitor how near the spacetime is to forming 
an apparent horizon, we employ the mass aspect func- 
tion, m, 



m(r,i) = -(l-l/a 2 ) 



(22) 



From Birkhoff's theorem, we know that an apparent 
horizon forms in the limit 2m /r — > 1. We note that 
Eqs. (|20j|21[) were used to compute the fluid equation 
source terms (|18H19p . and that flat space equations are 
obtained by setting 8 = and X = 1. 

With the conservation equations ([9^- (fT0f . the equa- 
tion of state (EOS) closes the system of hydrodynamic 
equations. Largely due to the extensive nature of our 
parameter space survey, we restrict the current study to 
continuum state equations (i.e. we do not use tabulated 
equations of state). The polytropic EOS 



is used only when calculating initial conditions for a star. 
Here, K is taken to be constant (isentropic condition) 
while r is the adiabatic index. After t = 0, we allow 
for the development of shocks and therefore only use the 
"ideal-gas" or "gamma-law" equation of state (|3|). Our 
initial neutron star models are approximated by solutions 
of the spherically-symmetric hydrostatic Einstein equa- 
tions, called the Tolman-Oppenheimer-Volkoff (TOV) so- 
lutions [H, [H, H(| • We simulate the stiffness of matter 
at super-nuclear densities by setting T = 2 in all of the 
calculations discussed here. As pointed out by Cook et 
al. [3l[ , the constant K can be thought of as the funda- 
mental length scale of the system, which one can use to 
scale any dynamical quantity with set values of (K, T) to 
a system with different values (K',T'). As with G and 
c, we set K = 1, thereby making our equations dimen- 
sionless. We note that for a specific value of T, the TOV 
solutions generically constitute a one-parameter family, 
where the central value of the fluid density, p c , serves as 
a convenient parameter. The ADM masses of stars for 
a TOV family governed by the ideal-gas EOS typically 
depend on the value of p c , with Madm — > for p c — > 0, 
and Madm achieving a global maximum at some value 
Pc — Pc- Stars with p c < p c are stable against radial per- 
turbations, while those with p c > p c are dynamically un- 
stable in radial perturbation theory. In the experiments 
described below, we generally work with stars that have 
central densities significantly less than p c . 

After the initial, star-like solution is calculated, an in- 
going velocity profile is added to drive the star to col- 
lapse. In order to do this, we follow the prescription 
used in [32j and [l9[. The method entails specifying the 
coordinate velocity 



_ dr u r 



(24) 



of the star. In general, the profile takes the algebraic 
form: 



Ug(x) = A Q (x 3 - B x) . 
The two profiles that were used in are 



U 2 (x) 



U, 



27 (7amp / 3 5x 
1 X 



10V5 



(25) 

(26) 
(27) 



P = Kp 



r 

o ) 



(23) 



where x = r/R+ and is the radius of the TOV solu- 
tion. Unless stated otherwise, the Ui profile will be used 
for all the results herein. The velocity profile is added 
consistently to the TOV solution by recalculating a and 
a via Eqs. (|20H2ip once the profile has been assigned to 
a given star. Further details concerning the calculation 
of initial data can be found in [l2], HH . 

As mentioned above, previous critical phenomena 
studies of perfect fluids have focused on models governed 
by the so-called "ultrarelativistic" EOS ©. This can 
be thought of as an ultrarelativistic limit of ([3]) , wherein 
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Po£ ^> Po or p ~ p Q e. In this limit, D becomes insignifi- 
cant and one is left with two equations of motion for the 
fluid, which can be easily derived from Eqs. lfT4 |) -(TT9 f by 
ignoring the EOM for D and setting E = r. The full 
expressions are given in 33]. In this paper, we only use 
the ultrarelativistic EOS in order to dynamically calcu- 
late ultrarelativistic Type II critical solutions. All other 
computations are performed using the ideal-gas EOS 



III. NUMERICAL TECHNIQUES & 
COMPUTATIONAL ISSUES 

In this section we discuss the numerical techniques we 
use to simulate the highly-relativistic flows encountered 
in the driven collapse of neutron stars. The simulations 
entail solution of a system of coupled, partial and or- 
dinary differential equations that describe how the fluid 
and gravitational field evolve in time. In section IIII Bl 
we also discuss an instability that generically appears 
in our calculations that are very close to the black hole 
threshold — this instability ultimately limits how closely 
we can tune any given family parameter to criticality. 

We use the Rapid Numerical Prototyping Language 
(RNPL) [U to handle check-pointing, input/output, 
and memory management for all our simulations. Sec- 
ondary routines are called to solve the fluid and geo- 
metric equations. We use second-order high-resolution 
shock-capturing (HRSC) methods to evolve the fluid. 
The discrete equations are derived using a finite volume 
approach and are detailed in Appendix [X] We generally 
use a Roe-type method (2(| HH as our approximate Rie- 
mann solver. However, particularly in the investigation 
of the instability mentioned above, we have sometimes 
used the Marquina flux formula [36| and Harten and Fly- 
man's [37j] entropy-fix for Roe's method, to compare with 
the basic Roe solver. For accurate and stable resolution 
of shocks, we use the minmod slope-limiter (38j to recon- 
struct the primitive variables at cell interfaces. We have 
also implemented the the linear MC [3!| and Superbee 
(40| limiters, as well as the high-order essentially non- 
oscillatory (ENO) scheme [4l|, but have found the min- 
mod limiter to provide the most stable evolutions near 
the threshold of black hole formation while still resolving 
shocks adequately. 

In order to track the continuously decreasing spatio- 
temporal scales typically seen in CSS phenomena, we 
use a nonuniform grid that refines as needed. The ori- 
gin, the point upon which the matter collapses, is the 
natural setting for the smallest dynamical scales and we 
therefore only refine the innermost region. Our imple- 
mentation was inspired by Neilsen [42| and is detailed in 
(33)]. The basic idea is to segment the discrete domain 
into three regions: an innermost uniform grid with the 
smallest grid spacing, Ar a , composed of N a cells, an adja- 
cent intermediate grid with cells of sizes Ar cx r, and 
an outermost uniform grid of N c cells. Refinement occurs 
when the maximum of 2m(r, t)/r is attained within N a /2 



cell widths from the origin. Interpolation is performed 
with a 3 rd -order ENO interpolation procedure written 
by Olabarrieta We usually set N a ~ 300 - 600, 

Nf, ~ 2N a , N c ~ 10 — 20 and adopt an initial value 
for Ar a such that the outer boundary lies at about 5-10 
times the initial radius of the star. 

Time integration is performed separately from the spa- 
tial discretization using the method of lines. Specifically, 
an explicit, two-step predictor-corrector technique, called 
Huen's method, is used to time-advance the ODEs that 
result from the spatial discretization of our time depen- 
dent PDEs. Discrete timesteps, At, for the ODE integra- 
tion are constrained to magnitudes given by At/Ar a < 
0.4, ensuring that the Courant-Friedrichs-Levy condition 
for our scheme is not violated. Additional details con- 
cerning the time integration are given in App. as well 
as in [33(]. Results from a shock tube test, which measures 
our code's ability to evolve discontinuities, and a conver- 
gence test involving the evolution of a self-gravitating 
distribution of fluid are presented in App. [U] 



A. Primitive Variable Calculation 

Since only the conserved variables are evolved by the 
HRSC schemes discussed above, the primitive variables 
must be derived from the conserved variables after each 
predictor and corrector step in order to compute fluxes f 
and source functions xp for the next evolution step. This 
involves inverting the three equations q = q(w) — given 
by the definitions of the conserved variables (JT5J) — for 
the three unknown primitive variables, w. While closed- 
form expressions for the inverted equations exist for the 
ideal-gas EOS, numerically solving the equations is far 
more efficient [4J|. At each grid point, we use a Newton- 
Raphson method to find the values of w that minimize 
the residuals of the conserved variable definitions (|15p . 
Instead of solving the full 3-by-3 system at each point, 
an identity function X — derived from Eq. (|15p — is used as 
a residual, making the solution process one-dimensional. 
This makes the procedure much more efficient. 

Our method for performing the inversion, which is dis- 
cussed further in App.lDl is based on one specially suited 
for spherically symmetric ultrarelativistic flows [251 X£l\ , 
and uses a residual function based on the definition of E: 

1(H) = HW 2 -t-D-P. (28) 

Here, H(r,t) = p a h is the enthalpy. In order to in- 
crease the accuracy of our computation of w, we use 
different methods for calculating the residual X and its 
derivative X' = dX/dH in different regimes (including 
both the ultrarelativistic and nonrelativistic limits). The 
"nonrelativistic" and "intermediate" methods originated 
from [25|, where flows in the ultrarelativistic limit 
were also studied. However, we have found that in the 
ultrarelativistic limit, where A = S/H — > 00, the inter- 
mediate method still gives imprecise results. This impre- 



cision can be traced to a loss of precision in the calcula- 
tion of the quantity 



1 - v (A) = 1 - v/l + l/4A 2 + 1/2 |A| 



(29) 



We thus expand all nonlinear expressions appearing in 
the conversion to primitive variables in powers of b = 
1/2 |A|, which yields results with increased floating-point 
accuracy. In the other limit, |A| <C 1, where the the 
system is nonrelativistic, we use expansions up to 0(A 9 ) 
that similarly reduce the influence of round-off errors. In 
practice, the ultrarelativistic regime is defined by an ad- 
justable parameter Anigh and the nonrelativistic regime 
by Alow For example, for all the results shown below, we 
used Afjigh = 10 2 and Al ow = 10 -4 ; these values ensure 
that the leading-order error terms in the ultrarelativis- 
tic and nonrelativistic expansions are below the intrinsic 
round-off error of the computations. 

Comparisons of the accuracy of our improved method 
and the work presented in [IE E3] are shown in Fig. [1] 
In order to estimate the relative error in the primitive 
variable calculation for each method, we seed each solver 
with guesses for w that are a fixed factor away from the 



(fc 



) - „,( fc ) 



^exLt (1 + 2 (fe) ), where 



true solution, 

is the exact solution, and denotes the k-th com- 

ponent of the state vector w. A constant set of seeds 
z (k) = {-0.2121,-0.0208941,-0.25971} are used in or- 
der to put all calculations on equal footing. Although 
the convergence of both methods does depend on the size 
of the our pre-specified initial values are generally 
further from the true solution values than they are in 
the context of an actual calculation. In addition, use 
of different seeds with magnitudes comparable to those 
given above yielded similar results, suggesting that these 
particular values of z^) are appropriately representative. 
For each method, once a best estimate for w is com- 
puted, the relative error for the solver is computed as 
(w - w oxact ) /wcxact- Fig. Q] shows the logarithm of the 
relative errors for the variables v and p (the trend of the 
error in P is similar to that of p ). The improved accu- 
racy of our new method is most noticeable in the compu- 
tation of v, as demonstrated by the first two columns in 
the figure. We note that our method can accurately cal- 
culate w for W > 10 3 and for all P, p tested, while the 
method described in [25| develops significant problems 
when W > 10 3 and P> Po . 

Even though the above methods improved the accu- 
racy of the primitive variable calculation, significant er- 
rors still remain for highly-relativistic flows (W > 10 5 ) 
where P and p Q differ by orders of magnitude — i.e. when 
P p or p 3> P. In these regimes, machine preci- 
sion limits the accuracy of the calculation since terms 
in X and X' become numerically insignificant (relative to 
other terms) even though their presence is essential to the 
computation of a solution. The effect from round-off er- 
ror in these regimes can easily be seen by plotting X{H) 
using different orders of numerical precision, which we 
have done using arbitrary precision arithmetic in Maple. 
In order to accurately calculate w in these regimes, one 
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FIG. 1: Comparison of the accuracy achieved with our prim- 
itive variable solver and that used in [25| ]. The first (third) 
column shows log 10 of the relative error between the exact 
value of v (£o) and the value obtained using the method de- 
scribed in [2g|, while the second (fourth) column shows log 10 
of the relative error between the exact value of v (p a ) and the 
value obtained from the method described in Sec. IIII Al In 
any given plot, the vertical (horizontal) axis is discretized into 
50 uniformly-spaced values of log 10 (P) (log 10 (p o )); the mini- 
mum (maximum) value used for both P and p is 2.5 x 10~ 17 
(10 6 ) . Each row contains plots of (P, p )-space calculations 
that were performed with a given value of W shown to the left 
of the row. The uniformly-spaced color map is shown at the 
bottom and indicates that the maximum (darkest shading) 
represents errors comparable to and exceeding 100%, while 
the minimum (lightest shading) represents errors comparable 
to and below machine precision. 



would need a new algorithm that performs better in these 
regimes, or use higher-precision arithmetic in the simula- 
tion. Fortunately, the improvements we have made seem 
to be sufficient for our current purposes. 



B. Instability at the Sonic Point in the CSS 
Regime 

In this section we provide a description of an instabil- 
ity that develops in our calculations in the vicinity of the 
sonic point in near-critical evolutions. When using the 
unmodified approximate Roe solver, this instability made 
it impossible for us to obtain consistent brackets about 
the critical parameter, p* , for \p— p*\ < 10~ 9 . This signif- 
icantly hindered our study, since we found that we needed 
to tune quite closely to the threshold solution in order to 
calculate an accurate value of the scaling exponent 7. 

An example of the instability seen in evolutions using 
primitive variable reconstruction is shown in Fig. [21 The 
conserved variable D is plotted in CSS coordinates T 



7 



and X defined later in the paper (Eqs. (j33|) and (j34|) . 
respectively). The last five frames show data from the 
last 5 time steps before the code crashes, while the first 
four frames are more distributed in T. Hence, we see that 
the feature at the sonic point exists for many discrete 
time steps before its growth produces a code crash. 
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FIG. 2: Conserved variable D(X,T) from the most nearly 
critical evolution obtained with the use of the approximate 
Roe solver without smoothing. X and T are defined by 
Eqs. (|34[) and (|33[1 . respectively. The dashed line indicates the 
location of the sonic point, X — 0. No refinement takes place 
during the period shown here, and Ar a ~ 1.55 x 1CP 7 . From 
left-to-right and top-to-bottom, the T values of the frames 
are -10.4109, -10.4977, -10.5916, -10.6938, -10.7822, 
-10.7823, -10.7824, -10.7825, -10.7826. The evolution 
started with a TOV star of central density p c = 0.05 that 
was perturbed using profile Ui (see (|26[0 with the overall am- 
plitude factor, ?7am P , tuned to produce near-critical evolution. 



The instability manifests itself in different ways, de- 
pending on the type of cell reconstruction used. For 
example when using the conserved variables to recon- 
struct the solution at the cell borders, we find that the 
conserved variables themselves remain smooth, but that 
each of the primitive variables exhibits persistent oscil- 
lations near the sonic point that span of order 2-4 grid 
cells. On the other hand, reconstructing with the primi- 
tive variables leads to smooth w but oscillations in q. A 
third and final reconstruction method was tried with the 
so-called characteristic variables, which are the advected 
quantities in the equations found by diagonalizing the 
quasi-linear form of Eq. (|14p . This method was signifi- 
cantly more diffusive but less stable than reconstruction 
with primitive variables. 

The instability is also sensitive to the slope lim- 
iter used, with Superbee and Monotonized Central- 



TABLE I: Asymptotic values of the fluid's characteristic 
speeds in the ultrarelativistic limit. The sonic point is lo- 
cated at r = r s . 



Characteristic Speed 


A(r < r s ) 


A(r > r s ) 


Ai 


< 1 


~ 1 


A+ 


~ 1 


~ 1 


A_ 


~ -1 


~ 1 



differenced (MC) limiters producing more spurious os- 
cillations than the more diffusive minmod limiter. We 
also found no improvements by varying the order of the 
ENO reconstruction from 0(Ar 2 ) through 0(Ar 10 ). 

In terms of ruling out potential sources of the problem, 
we have ensured that the regridding procedure is not re- 
sponsible for the instability. To accomplish this, we first 
evolved a system that was tuned near the critical solu- 
tion. We extracted the grid functions at a specific time, 
t, before the appearance of instability, and interpolated 
them onto a new grid fine enough so that no further re- 
finement would be required in the subsequent evolution. 
The data was allowed to evolve from this time, and the 
instability developed in the same manner and at the same 
time, t, as in the original run. 

Moreover, we find that the instability does not "con- 
verge away." We tuned the initial data towards criticality 
for three different levels of refinement, where refinement 
was done locally so that Ar; (r) = 2 Ar/+i (r) for all r, and 
I is the "level" of refinement. We find that as I increases 
the oscillations associated with the instability do not sig- 
nificantly change in magnitude and remain confined to 
approximately the same number of grid cells. Also, the 
solutions eventually diverge at the sonic point in all cases. 

In order to describe the likely source of the instability, 
we first need to provide a better description of the near- 
critical solution. When the initial data has been tuned 
close to the critical solution at the threshold of black 
hole formation, the behavior near the origin is self-similar 
up to the sonic point, r s , where the flow velocity equals 
the speed of sound, c s (|A12|) . For these solutions the 
fluid becomes ultrarelativistic — e.g. P » p — for r < 
r s and we expect that c s (r) — ► 1 in that region. Also, 
from previous ultrarelativistic studies using r = 2 such as 
[E0; we expect that v — ► 1 for r > r s . Thus, about the 
sonic point, the characteristic speeds (IA6[) should take 
the values given in Table |TJ 

In fact, this is exactly what we find when using the 
ideal-gas state equation, as seen in Fig. [3] and Fig. 0J In 
Fig. [4] we see that P p a within the self-similar region, 
but that P(r) < p {r) for r > r s . 

From these plots we also find that the transition from 
the ultrarelativistic regime to the exterior solution — 
defined by r > r s , and characterized by an absence of 
self-similarity and decay to asymptotic flatness — is quite 
abrupt. For instance, the discontinuity in A_ is resolved 
by only a few grid points, signifying the presence of a 
shock which can also be seen for r ~ r s in the plots of 
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FIG. 3: Characteristic speeds of the fluid for the most nearly 
critical solution obtained with the approximate Roe solver 
without smoothing. The wave speeds are plotted here as 
functions of the self-similar coordinate X, and are shown at 
T = —10.6938. A closer view of the characteristic speeds 
near the sonic point is shown as an inset in the lower-right of 
the plot, revealing the severity of the discontinuity in A_ as 
discussed in the text. 



P(r) and p a (r) shown in Fig. 2] Since A_(r < r a ) < 
and A_(r > r s ) > 0, the discontinuity represents a point 
of transonic rarefaction. Also, the shock appears to be 
an expansion shock, which is entropy-violating, since it 
travels into a region of higher pressure and density. The 
reason why we can have an entropy-violating shock de- 
velop is because it is coincident with a change in curva- 
ture: as the high-pressure matter leaves the confines of 
the potential, it freely expands and its internal energy is 
converted into bulk kinetic energy. 

LeVeque states in [45] that the Roe solver can lead 
to the wrong Riemann solution at transonic rarefactions 
(in flat spacetime) since the linearization that the Roe 
solver performs on the EOM leads to a Riemann solution 
with only discontinuities and no rarefaction waves. He 
illustrates this point in [46| using a boosted shock tube 
test that makes the rarefaction transonic. Other failures 
of Roe's method that are attributed to its linearization 
have been shown by Quirk [47[ , and by Donat et al. [48| 
where an unphysical "carbuncle" forms in front of a rel- 
ativistic, supersonic jet. 

To see whether Roe's method contributes to the 
instability, we implemented the Marquina method 
and the entropy-fix for the Roe scheme due to 
Harten and Hyman [371 ] . A comparison between 
the three methods is shown in Fig. where we 
have evolved a shock tube problem that emulates the 
fluid state about the sonic point of near-critical so- 



FIG. 4: Pressure and rest-mass density of the most nearly 
critical solution obtained with the approximate Roe solver 
without smoothing. Both quantities are plotted as a func- 
tion of the self-similar coordinate X, and are shown at T = 
— 10.6938. Interior to the sonic point, X = 0, the fluid is 
clearly in the ultrarelativistic limit, with P/p a ~ 10 4 . How- 
ever, beyond the sonic point, the flow is not ultrarelativistic — 
in fact P < p in most of the domain exterior to X — 0. A 
closer view of the distributions near the sonic point is shown 
in the inset, and more clearly illustrates the formation of an 
expansion shock as discussed in the text. 



lutions. The initial conditions used for this test 
are {p L ,v L ,P L } = {l.O x 10 3 , -0.3, 1.0 x 10 6 } and 
{pR,v R ,P R } = {0.3,0.9994, 1.0}. These values are such 



that, initially, A +L ~ 0.9995, A H 



0.99998, A_ 



—0.99987, and A_ H ~ 0.98296, closely approximating 
the values listed in Table |U The Roe, Marquina and 
Roe with entropy-fix evolutions each used 400 points in 
the entire grid of domain < x < 1 (only a portion of 
the grid is shown here), and both used a Courant factor 
of 0.4. The exact solution was obtained from the Rie- 
mann solver provided in |49[ with 1000 points, using the 
same range in x and same initial conditions as the Roe 
and Marquina runs. The Marquina method and entropy- 
fixed Roe method produce similar — yet more diffused — 
solutions than the exact solver, while the basic Roe solver 
severely diverges from the exact solution near the tran- 
sonic rarefaction during the first few time steps. Even 
though the initial divergence of the Roc solution eventu- 
ally vanishes, a relic feature still exists and propagates 
away from the center. If we were to reverse the evolution 
of the Roe solution shown here, the sequence would be 
reminiscent of how the instability in D grows near the 
sonic point of near-critical solutions (Fig. [2]) . 

This shock tube test suggests that our use of Roe's 
method may underlie the observed instability. In order 
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FIG. 5: One-dimensional, slab-symmetric shock tube test to 
simulate the discontinuity observed near the sonic point of 
near-critical solutions. The rest-mass density p (x,t) com- 
puted using different Riemann solvers is plotted versus x in 
each constants frame. Solution time is shown in the upper- 
right part of each frame. The solid line without points is 
an exact Riemann solution, the dashed line with triangles 
corresponds to the solution obtained with the approximate 
Roe solver, squares represent the solution from Marquina's 
method, and circles are from the Roe scheme with entropy 
fix. 



to address this possibility, we implemented both the Mar- 
quina solver and the entropy-fix in the general relativistic 
code and tuned towards the critical solution. We were 
able to tune to \p — p*\ « 5 x 10~ n with Marquina's 
method, which is approximately a factor 10 2 closer to p* 
than we reached with Roe's method. Also, the "bump" 
at r s developed at a later T with Marquina's method. 
However, the use of Marquina's flux formula did not 
completely solve the problem since evolutions using it 
also eventually succumb to the instability, preventing 
us from tuning beyond \p — p*\ 5 x 10 -11 . More- 
over, Roe's method with the entropy-fix, performed more 
poorly than when used without the entropy-fix, impeding 
any further tuning past \p — p* \ ~ 2 x 10~ 8 . 

A possible explanation for the failure of the entropy- 
fixing methods to stabilize near-critical evolutions may 
involve the fact that they do provide an entropy-fix. As 
mentioned before, these methods would evolve the criti- 
cal solution's expansion shock to a rarefaction fan if no 
source was present (see Fig. [5J. However, the source in 
the EOM and/or the geometric factor in the flux must 
be what keeps the expansion shock from being dissipated 
away since the shock is a part of the critical solution. 
The competition between the approximate Riemann so- 
lutions of the entropy-fixing methods and the spacetime 



may exacerbate the instability. However, it is unclear 
why Marquina's method is more successful than Roe's 
method with an entropy-fix. 

Even though Marquina's method provides a marked 
improvement over Roe's method, it does not eliminate 
the sonic-point instability. A first attempt to explicitly 
dissipate the instability involves applying artificial vis- 
cosity in the region. We follow Wilson's [50j artificial 
viscosity method and set P — > P + Q in f alone, where 



(30) 



and c av is a user-specified parameter. However, the insta- 
bility worsens as v — > 1, and, since Q is not proportional 
to W like other terms in f , Q becomes irrelevant as the 
flow becomes more relativistic. As a result, moderate val- 
ues of c av lead to insignificant changes in behavior near 
the sonic point, while extremely large values tend to am- 
plify the instability and induce additional high-frequency 
modes near r s . 

Since we find that the instability near r s becomes more 
severe as A_ became more discontinuous, our second at- 
tempt to control the blow-up entails smoothing the con- 
served variables about r s at every predictor/corrector 
step of the fluid update. The smoothing is done once 
the A_ discontinuity develops. Specifically, we start to 
use the smoothing procedure when \p— p*\ < 10~ 8 — 10~ 9 , 
and at times when A_ begins to be resolved over approx- 
imately 10 or fewer zones. We can use the same time, 
t s , to begin smoothing for all runs since the evolution for 
t < t s is almost identical for all near-critical values of p. 
The smoothing is performed over the first contiguous set 
of points, r™, that satisfy -A min < A_(r| m ) < A min , 
where A™ ln is some adjustable parameter which we set to 
0.95. The smoothing operation replaces the quantity qj 
with (q j+ i + <?j-i)/2. 

We also find that the instability worsens as the num- 
ber of points between the origin and the sonic point de- 
creases, as occurs in those cases where the solution dis- 
perses from the origin instead of forming a black hole. 
Our ability to follow evolutions through to dispersal is 
necessary for calculation of the scaling exponent, 7, since, 
as described in the next section, one way of estimating 
7 involves measuring the scaling of the global maximum, 
r max , of the stress energy trace, T(r, t) = T a a , as a func- 
tion of \p — p*\ . We find that T max is generally attained 
at a time when the fluid is beginning to disperse. Conse- 
quently, we found it necessary to refine the grid whenever 
the discontinuity or max,. (2m /r) reaches r ~ r a j1. 

The diffusion introduced by the smoothing allows us 
to further tune toward the critical solution, eventually 
to \p — p*\ ~ 5 x 10~ 12 . However, we are still unable to 
calculate the global maximum of T, T max , for the most 
nearly critical runs even though we can identify them as 
being dispersal cases. The minimum value of \p — p*\ 
for which we can calculate T max is about 5 x 10~ 10 , as 
illustrated in Fig.[5J This is far smaller, however, than we 
can achieve without smoothing or with any other method 
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TABLE II: Star solutions in which we observed Type II be- 
havior, and the minimum black hole masses we were able to 
form from them. We denote the mass of the smallest black 
hole found for a given p c by min(MBH), M* = M ir (p c ) is the 
mass of the initial star solution, and min \p — p* /p is the rel- 
ative precision reached in p* per star. The final columns lists 
the critical parameter values we obtained. 



central density p c = 0.05. 



p c min(A/ B H)/M* 


min \p — p*\ lp 


P* 


0.01 1 x 10~ b 


2 x 10 _y 


0.889 


0.02 6 x 10~ 7 


1 x 10~ 9 


0.746 


0.03 3 x 10~ 7 


5 x 10~ 10 


0.634 


0.04 6 x 10~ 8 


2 x 10 -11 


0.543 


0.05 2 x 10~ 8 


6 x 10" 12 


0.469 


we have tried. Surprisingly, smoothing q about the sonic 
point did not make the Marquina evolutions any more 
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IV. SIMULATION RESULTS 
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In this section, we study the Type II, CSS critical solu- 
tion found at the black hole- forming threshold of the pa- 
rameter space described in If not otherwise stated, 
the results in this section use U\ (127p for the initial ve- 
locity profile, and the overall amplitude, £/ amp of the 
velocity profile for the tuning parameter, p. As previ- 
ously mentioned, the stars that we able to drive to a 
Type II black hole threshold generally have central den- 
sities, p c , significantly smaller than the maximum value, 
p c , along the stable branch, which for T = 2 is 0.32 in 
our units. Although we were generally able to form black 
holes from stars with an initial rest-mass central density 
greater than p" lm = 0.007, we have closely tuned to- 
wards critical solutions for only a handful of such initial 
states. (We were unable to form black holes from stars 
with p c < p™ 111 using an initially-ingoing velocity profile.) 
In Table [D] we list the central densities of the stars for 
which Type II behavior was actually observed, and quan- 
tify how close to the critical value we were able to tune. 
The instability described in Sec. IIII 51 limited the tuning 
in all instances. 

From Table HTl it is clear that the instability's effect on 
our ability to find the critical parameter increases with 
decreasing p c . This is most likely due to the fact that 
sparser stars require greater in-going velocities in order 
to collapse, giving rise to more relativistic and, conse- 
quently, less stable evolutions. We note, however, that 
our results represent great improvement over the preci- 
sion obtained in [lj|; the smallest black hole attained in 
that study was min(MBH")/M* ~ 10~ 2 . The success of 
our code is most likely due to our use of adaptive/ variable 
mesh procedures and the great lengths we went to com- 
bat the sonic point instability. 

Unless otherwise stated, and for the remainder of the 
section, we focus on behavior seen with the star having 



FIG. 6: Scaling behavior for supercritical — or black hole 
forming — solutions. The top plot illustrates how the points 
from a series of supercritical runs follow the scaling law for 
the black hole mass (pTJ , while the bottom plot shows how the 
data deviate from our best fit to this scaling law. The two 
dotted lines delineate the data used in making the best fit; 
this data is plotted separately in Fig. [7] Black holes were as- 
sumed to have formed when max r (2m(r,t)/r) > 0.995. The 
gaps between some of the points represent those runs that 
crashed before max r (2m/ r) reached this value. Smoothing 
was used for In \p — p*\ < —19.3, which is also where we start 
our fit. These runs used p c = 0.05, U = Ui and an initial grid 
defined by {N a ,N b , N c , Ar a , level} = {300, 500, 20, 0.005, 0}. 

To demonstrate the scaling behavior of Mbh , we show 
in Fig. [5] In (Mbh) versus In \p — p*\ for a wide range of 
supercritical solutions. The slope of the trend is equal to 
the scaling exponent, 7. From the figure, we can clearly 
see that the scaling law provides a good fit only in the 
limit p — > p* as expected [liij . The jump seen at In \p — 
p*\ ~ —10 represents the point at which the fluid enters 
a dynamical phase where the center part of the star has 
enough kinetic energy to dominate the effective potential 
energy, whose magnitude is set by p . In this regime, the 
fluid then follows a CSS-type evolution. 

In addition, Fig. [5] is meant to illustrate problems in 
our calculations associated with the coordinate singu- 
larity that inevitably develops in our Schwarzschild-like 
coordinate system in super-critical evolutions. Compu- 
tations were run for a set of parameter values, pi, dis- 
tributed uniformly in In \p — p*\ — any gaps in the plotted 
data thus represent instances where our code crashed pre- 
maturely. We note that the flow velocity becomes discon- 
tinuous and nearly luminal when black hole formation is 
imminent, and this seems to amplify the instability men- 
tioned in Sec. IIII Bl This results in the evolution halting 
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ln(p - p*) 

FIG. 7: Best-fit for the scaling behavior of black hole masses 
near the threshold. The top plot shows calculated masses and 
the fitting line, while the bottom plot shows the deviation 
between the two. The scaling exponent for this fit, which is 
simply the slope of the line, is 7 = 0.94. 

before max r (2m(r, t)/r) exceeds its nominal threshold of 
0.995. However, for a set of parameter values delineated 
in the figure by dashed lines, we were able to find a good 
fit to a scaling law. For that subset of data, the fit, as 
well as the data's deviation from the fit, are shown in 
Fig. [71 One measure of how well the black hole masses 
are described by a relation of the form (TTJ) is that devia- 
tions from the fit are small and apparently random. The 
slope of the trend yields an estimated scaling exponent 
of 7 = 0.938. 

As mentioned in the previous section, to obtain an- 
other estimate of the scaling exponent, we calculate how 
the global maximum, T max , of the stress-energy trace 
scales as p — > p*, using subcritical computations. As 
Garfinkle and Duncan [51] pointed out for the case of 
spherically-symmetric massless scalar collapse, the global 
maximum of the Ricci scalar should be proportional to 
the inverse square of the fundamental length scale of the 
self-similar solution. Hence T max for near critical solu- 
tions below the threshold should follow the scaling law: 

T max c< b-pT 27 ■ (31) 

Using T instead of the Ricci scalar is computationally 
more expedient since it does not require the calculation 
of second-order space and time derivatives of the metric 
functions. 

By determining 7 from a plot of T max versus \p — p*\, 
in addition to a fit to Eq. (TT|), we can get an estimate 
of the systematic errors in our estimation of 7 for both 
methods. Estimation of 7 from the scaling of T max also 
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FIG. 8: Scaling behavior in T max for subcritical solutions, 
i.e. those not forming black holes. The line shown here is 
the best-fit for the expected scaling law (|3ip . using data 
only from the solutions closest to criticality. These runs 
used p c = 0.05, U — Ui and an initial grid defined by 
{N a , N b , No, Ar a , level} = {300, 500, 20, 0.005, 0}. 

has the advantage that, in the limit p — > p*, the code 
is more stable for subcritical, rather than supercritical, 
evolutions. The scaling behavior for T max can be seen 
in Fig. [5] where lnT max is plotted versus In \p — p*\. The 
solutions far from criticality seem to smoothly asymptote 
toward the critical regime. The line shown in this plot 
uses only those points in the regime that provide the 
best linear fit; a closer view of the points used in the fit 
are shown, for instance, in Fig. [T2] Since the slope of 
the line now represents —27 (|31[) . we find from this fit 
that 7 — 0.94, which agrees with the value found from 
the scaling of Mbh(p) to within the estimated systematic 
error in our computations. 

Although our calculated scaling exponents match well 
to results previously obtained for the ultrarelativistic 
fluid with r = 2, this does not necessarily say how well 
the ideal-gas critical solutions compare to the ultrarel- 
ativistic ones in detail. To obtain the ultrarelativistic 
critical solutions, we let an adjustable distribution of ul- 
trarelativistic fluid free-fall and implode at the origin; 
specifically, the initial data for the fluid is set so that 
r(r, 0) is a Gaussian distribution and S(r, 0) = 0, and 
the amplitude of the Gaussian is used as the tuning pa- 
rameter. The scale-free functions from the critical so- 
lutions of the velocity-induced neutron star system and 
the ultrarelativistic system are shown in Fig. [5J Here, 
the quantity to = uj(r,t) is another scale- free function 
determined from metric and fluid quantities via 

u) = 4irr 2 a 2 p. (32) 
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FIG. 9: Three scale-free quantities of near-critical solutions in 
self-similar coordinates for the ideal-gas system (dashed line) 
and the ultrarelativistic system (solid line). 

In order to make the comparison between the two so- 
lutions, the grid functions were transformed into self- 
similar coordinates T and X: 

T = In (T * -To) , (33) 




where Tq is the elapsed central proper time, 




Tq is the accumulation time of a given critical solution, 
and r e (t) is the location of the sonic point. Since we find 
that computing a smooth r s (t) near the critical point is 
difficult, we typically use X a , 

X a = In (y—^j > (36) 

as our self-similar radial coordinate. Here, ?" 0max (i) is the 
position of the local maximum of a(r,i) that lies closest 
to r = 0. 

Our results indicate that the ideal-gas system does 
asymptote to the ultrarelativistic self-similar solution in 
the critical limit. While the ultrarelativistic fluid enters a 
self-similar phase shortly after the initial time, the ideal- 
gas solution generally tends toward the critical solution 
relatively slowly, then eventually diverges from it. The 
agreement between the ideal-gas and ultrarelativistic so- 
lutions improves as p — > p* , as expected, and Fig.[9]shows 



FIG. 10: Deviation over time of those quantities displayed in 
Fig- EI Here, ll/H denotes the ^2-norm of the function /. The 
deviations for a (solid line), u> (dotted line), and v (dashed 
line) are shown. The fe-norms of these differences are com- 
puted at every time satisfying X a < 2, and then logarithms 
of those norms are plotted as a function of self-similar time 
T. Note that the sense of physical time is opposite to that of 
T; that is, T — » — oo as the solution approaches the accumu- 
lation time. As the evolution proceeds from the initial time, 
the two solutions asymptote toward each other. For T < —13, 
the deviation between the two solutions increases as the ideal- 
gas near-critical solution departs from the asymptotic critical 
solution and eventually disperses from the origin. 



profiles at a time when the difference between the solu- 
tions was minimized. The -^-norms [55j of the deviations 
between the ideal-gas and ultrarelativistic scale- free func- 
tions plotted over time are shown in Fig. 1101 it can be 
easily gleaned from this figure that the minimum of the 
average deviations occurs at approximately T = —13.1, 
which is the time at which we have displayed the profiles 
in Fig. [5] Also, Fig. [TD] graphically illustrates how the 
ideal-gas solution asymptotes exponentially — in central 
proper time, T — to the ultrarelativistic critical solution 
at early times. The deviations for the three functions 
seem to have the same qualitative trend, indicating that 
metric and fluid quantities asymptote to their ultrarela- 
tivistic counterparts. 

This exponential approach of the ideal-gas solution to 
the self-similar solution is better seen in the comparison 
of time sequences of ui extracted from the ideal-gas and 
ultrarelativistic computations, as shown in Fig.[TT] Here, 
^uitra has already attained a self-similar form at the be- 
ginning of the displayed sequence, while Widcai becomes 
self-similar at later times, and remains self-similar only 
for a time interval AT ~ 6. 
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FIG. 11: Time sequences of ui for the most nearly critical 
solutions obtained with the ideal-gas EOS (dashed line) and 
the ultrarelativistic EOS (solid line). Both functions have 
been transformed into self-similar coordinates, based upon 
their respective accumulation times and respective values of 
r a m ax- Approximate values of T are shown in the upper- 
left corners of the frames. Note that the ultrarelativistic uj 
is varying slightly frame-to- frame, contrary to appearances. 
Relative to the intrinsically ultrarelativistic solution, it takes 
more time for the ideal-gas solution to become self-similar 
since the length scale set by p a in the latter case only becomes 
insignificant for P/p a 2> 1. 



A. Universality and Consistency 

This section describes several numerical experiments 
primarily designed to ensure that the results presented 
above are not artifacts of the computational techniques 
used. These computations also provide a measure of the 
systematic error in our calculation of 7. Moreover, in or- 
der to check previous claims that critical solutions gen- 
erated from perfect fluid configurations having the same 
adiabatic index T may not reside in the same universality 
class, we also measure 7 for different initial conditions, 
while keeping T constant. When making any compar- 
isons, the methods, parameters, and initial data used 
to produce Figs. [SH5] will be referred to as the "origi- 
nal" configuration. A tabulation of the values of 7 and 
p* calculated from the different simulation configurations 
discussed in this section is given in Table IIIII 

The effect on the scaling behavior due to the floor 
value, 5, used for the fluid (see App. is estimated 
first. Since the magnitude of the floor is set in an ad hoc 
fashion — i.e. without any physical basis — it is crucial to 
verify that any results are independent of it. To test this, 
we replicated the original runs using different values of 
the floor while keeping all other parameters fixed. The 



FIG. 12: Scaling behavior in T max near the critical solution 
for runs using different values of S. The scaling behaviors 
are shown for the original configuration (circles,solid line), a 
configuration with 10 2 times the original floor value (squares, 
dotted line), and one with 10 4 times the original floor value 
(triangles, dashed line). The scaling exponents 7 for these 
runs are listed in Table ITTT1 



scaling behavior obtained using different floor values is 
illustrated in Fig. [121 The dotted and dashed lines cor- 
respond to floor values that are factors of 10 2 and 10 4 , 
respectively, larger than the original configuration, which 
itself used S = 2.5 x 10 -19 . The minimal influence of the 
floor on solutions in the critical regime is clearly seen by 
the fact that all points follow nearly the same best-fit 
line. In fact, Table [TTTI indicates that all estimated val- 
ues of 7 agree to within ~ 0.5% and that all estimates 
of p* coincide to within 0.0005%. The deviations of the 
calculated sets {In (T max ) , In \p — p*\} from their respec- 
tive best-fit lines for the different floor values even follow 
the same functional form, suggesting that the observed 
"periodic" deviations from linearity are not due to the 
floor. 

The absence of any dependence on the floor is not too 
surprising since the component of the fluid that under- 
goes self-similar collapse is never rarefied enough to trig- 
ger use of the floor. For instance, at a time when the cen- 
tral part of the star begins to resemble an ultrarelativistic 
critical solution, the minimum values of {D, II, $} within 
the sonic point are, respectively, {~ 10 2 , ~ 10 3 , ~ 10 3 } — 
far above the typical floor values used. Only for r > 
is the floor activated, and dynamics in this region can- 
not affect the interior solution once self-similar collapse 
begins due to the characteristic structure of near-critical 
solutions (see Table fl}. 

The effect from the Riemann solver used on the scaling 
behavior is seen in Fig. [T3] We find that the scaling be- 
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FIG. 13: Comparison of the scaling behavior in T max obtained 
with two different Riemann solvers. The "Smoothed Roe" line 
corresponds to the original calculations. The other (dotted) 
line was generated using the Marquina method, with other 
computational methods and parameters identical to the orig- 
inal calculations. The scaling exponents, 7, for these runs are 
listed in Table Hill 



havior of T max from the two methods is remarkably close. 
Even though the Roe method with smoothing allows us 
to determine In (T max ) for smaller values of In \p— p*\, the 
deviations from the best-fit of the two data sets are of the 
same order of magnitude for common values of In \p — p* \ . 
From Tabic IIII1 we see that the respective values of 7 
agree to within 0.3% and that values of p* agree to within 
0.001%. These differences are quite small — comparable 
to those found as a result of varying the floor. Hence, we 
conclude that the choice in Riemann solvers has little, if 
any, effect on the computed scaling behavior, indicating 
that the smoothed approximate Roe solver is adequate 
for our purposes. 

When using finite difference methods, it is vital to ver- 
ify that the order of the solution error is the same as the 
order to which the derivatives are approximated by differ- 
ence operators. For example, our HRSC scheme should 
be 0(Ar 2 ) accurate in smooth regions and O(Ar) near 
shocks, so we should expect this scaling behavior of the 
error as Ar is varied. First, we wish to see if our esti- 
mate for 7 converges as the grid is refined. Figure 1141 
shows a plot of In (T max ) versus In \p — p*\ for the origi- 
nal configuration, along with others computed at higher 
resolutions. We first see that the three distributions fol- 
low lines of approximately the same slope (and which 
are thus shifted vertically relative to one another by con- 
stant amounts) while the deviation of the best-fits seems 
to increase slightly with resolution. Also, we can see 
that an increase in resolution permits us to follow the 



FIG. 14: Scaling behavior in T max near the critical solution 
for runs using different levels of resolution. The runs were 
made with p c — 0.05, U = Ui, and the solid line with cir- 
cles was generated from runs using the original configuration. 
The level = 1 (squares, dotted line) and level = 2 (triangles, 
dashed line) runs, respectively, used computational grids that 
were locally 2 times and 4 times as refined. The scaling ex- 
ponents, 7, for these runs are listed in Table ITTTI 



collapse through to dispersal for solutions closer to the 
critical threshold, allowing for the scaling law to be sam- 
pled at smaller In \p — p*\ . Even though the deviations 
from the best-fits for I = 1, 2 are quite small compared 
to the typical size of ln(T max ), it is a little worrisome 
that they are larger than those from the lowest resolu- 
tion runs. However, this behavior can likely be attributed 
to the sonic point instability and the smoothing proce- 
dure used to dampen it. In particular, the "bump" at r s 
sharpens with increasing resolution spanning a roughly 
constant number of grid cells (see Sec. MI 51 for more de- 
tails). Consequently, the impact of the instability on the 
solution may also increase with decreasing Ar, since the 
discretized difference operators will generate increasingly 
large estimates for spatial derivatives in the vicinity of the 
sonic point. In addition, the smoothing operation is al- 
ways performed using nearest-neighbors, so the smooth- 
ing radius physically shrinks with resolution, diminishing 
the impact of the smoothing. 

In order to verify that the code is converging in the 
self-similar regime, we computed independent residuals 
(i.e. applied discretizations distinct from those used in 
the scheme used to compute the solution) of the Hamil- 
tonian constraint (f20)) and slicing condition (]21[) for the 
three levels of resolution discussed previously (Fig. 1T5)) . 
The overlap of the scaled residuals seen in the figure in- 
dicates 0(Ar 2 ) convergence. Note that the smoothing 
procedure has not been used to calculate the solutions 
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FIG. 15: Logarithm of scaled, independent residuals of the 
Hamiltonian constraint (|20[) an d slicing condition (|21[) for 
three levels of resolutions calculated from solutions in the 
self-similar regime. The dotted (dashed) lines are from a run 
which used 2 (4) times the local spatial and temporal reso- 
lutions of the original run, which is represented by the solid 
lines; the dotted (dashed) residual was scaled by a factor of 4 
(16) in order to make the 0(Ar 2 ) convergence of the solution 
more apparent. Each distribution is from a solution that has 
been tuned to In \p — p*\ ~ —19 with respect to the value of 
p* for each resolution, and function values at every tenth grid 
point are shown. The physical velocity of the fluid for the 
I — run is shown in the bottom frame in order to facilitate 
comparison of features in the independent residuals to those 
in the solution. 

shown here. We see that the scaled residuals have simi- 
lar magnitudes in all regions except those that have been 
processed by shocks, namely X a = [0, 4.5],~ 7.8, ~ 9.4. 
Because the self-similar solutions are converging at the 
expected rate, we surmise that the variations observed in 
7 for the three resolutions does not indicate a problem 
with convergence, but demonstrates the effect of trunca- 
tion error and/or the smoothing procedure on the scaling 
behavior. With only three levels of resolution, it is hard 
to make definite claims as to whether 7 is or is not con- 
verging to a particular value. Even so, the standard devi- 
ation of 7 determined from the three evolutions is about 
1% of their mean, suggesting that the variation is not 
significant. In fact, it is comparable to the 2% standard 
deviation found from the simpler ultrarelativistic perfect 
fluid studies of [7[. 

The final comparison entails varying the physical ini- 
tial conditions of the system to investigate the universal- 
ity of the critical phenomena computed with the ideal-gas 
EOS. The primary constituents of our model are the ini- 
tial star solution and the form of the perturbation with 
which we drive the star to collapse. Hence, we choose 



FIG. 16: Scaling behavior in T max for several families of initial 
data. The "Original" date (solid line) is as before, the dotted 
line with squares shows the scaling behavior for runs that used 
a different initial velocity profile, U = U%, and the dashed 
line with triangles was made from runs with a different TOV 
solution, with p c — 0.0531. The scaling exponents, 7, for 
these runs are listed in Table iHll 

to perform sets of runs to measure the scaling law us- 
ing 1) a different initial star solution and 2) a different 
functional form of the initial velocity profile. The scaling 
behaviors of In (T max ) versus In \p — p*\ for these differ- 
ent configurations are compared to the results from the 
original configuration in Fig. 1161 For the data computed 
using a star of central density p c — 0.0531, the only dif- 
ferent aspect is the initial star solution. This particular 
value of p c is chosen since it is near the transition from 
Type II to Type I phenomena discussed in detail in [l2| . 
The second initial data set uses U2 (|27p for the initial 
profile of the coordinate velocity. Naturally, we see that 
the three data sets are vertically shifted relative to one 
another since each set evolved from significantly different 
profiles of mass-energy — the details of the initial data set 
the scale for T max for specific values of In \p— p*. However, 
only the slopes of the curves are relevant for estimating 
7- 

From the values listed in Table HTT1 we see that 7 varies 
more significantly with the particular star solution used 
than with the form of the velocity profile. In fact, we are 
able to tune closer to the critical solution with the more 
compact star, a trend that can also be seen in Table [Til 
Nonetheless, the scaling exponents computed using the 
two distinct initial star configurations agree to within 2%. 

The change in the initial velocity profile only affects the 
computed value of 7 by 0.04%. This suggests that other 
methods of perturbation would also yield close to the 
same value. The concordance of results from these three 
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TABLE III: Scaling exponents 7 and critical parameters p* 
computed from fits to the expected scaling behavior in T max . 
Scaling exponents were obtained from runs using initial star 
solutions with different central densities (p c ) and using differ- 
ent floor magnitudes (5), levels of refinement (I), and velocity 
profiles (U). The runs labeled "Roe" use the approximate Roe 
solver with smoothing, the "Marquina" run used the Mar- 
quina flux formula, and the "Ultra-rel." scaling exponent was 
computed from our results involving the collapse of Gaussian 
profiles of ultrarelativistic fluid. 
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different families of initial data imply that universality of 
critical solutions is maintained for perfect fluids governed 
by an ideal-gas EOS, at least for the case T = 2. It would 
be interesting to see whether these properties hold with 
even more realistic equations of state, as well as for other 
values of T. 



B. Final Determination of 7 

Using the calculated values of 7 from the various meth- 
ods, floor sizes, and grid resolutions, we are able to pro- 
vide an estimate of the systematic error inherent in our 
numerical model. Further, by assuming that the univer- 
sality is strictly true, we can even use the variation in 7 
computed from the different initial data families for this 
estimate. Taking the average and calculating the stan- 
dard deviation from all of the values for the ideal-gas 
EOS listed in Table Hill we estimate a scaling exponent 
value of 

7 = 0.94 ±0.01. (37) 

This is in agreement with the value of 7 computed from 
the black hole mass scaling fit shown in Fig. [7] 

In addition, we can compare our final estimate of 7 
to values previously found for the ultrarelativistic fluid. 
As already mentioned, 7 was measured at three different 
refinement levels in [7[, and a value 

7< 0.96. (38) 

was quoted. 

Instead of solving the full set of PDEs, 7 can also be 
found by solving the eigenvalue problem that results from 
performing first order perturbation theory about the CSS 



solution. This was done in two ways in [f| : using the com- 
mon shooting method, and solving the eigenvalue prob- 
lem directly after differencing the equations to second or- 
der. The scaling exponents calculated were, respectively, 
7 = 0.9386 ± 0.0005 and 7 = 0.95 ± 0.01. 



C. Possible Presence of a Kink Instability 

As mentioned in Sec. IIII B[ we witness an instabil- 
ity near the sonic point of solutions tuned close to the 
threshold. After thorough numerical experimentation, 
we are still left uncertain about its genesis. One possi- 
bility is that it has physical origin. For instance, Harada 
[52| reported the presence of an unstable kink mode 
for spherically symmetric ultrarelativistic perfect fluids 
when r > 1.889. The unstable kink mode manifests it- 
self as an ever-steepening discontinuity in p at the sonic 
point, r s , that diverges in finite proper time. A mild, seed 
discontinuity at r s is necessary for the mode to grow, 
but — since it diverges in finite time — minute disconti- 
nuities inevitable in discretized solutions, for instance, 
are expected to be sufficient to manifest the instability. 
Physically, the "seed" discontinuity could be the scale at 
which the continuum approximation of hydrodynamics 
breaks down, i.e. at the particle scale. 

In our nonlinear PDE solutions that use the ideal-gas 
EOS, we do in fact find a growing discontinuity in P and 
p a at r s (Fig. 2]), and this is precisely where our instability 
develops. 

However, as was the case in Q, our solutions of the 
PDEs for an explicitly ultrarelativistic fluid do not de- 
velop instabilities in the vicinity of r s for T = 2. One 
possible reason for this could be the fact that in this case 
the flow is never transonic, i.e. there is no sonic point 
since c s — 1, and the flow can never attain this velocity. 
For the ideal-gas fluid, there is a sonic point, as c s < 1 
(albeit arbitrarily close to 1) since p /P will always be 
non-zero during a numerical evolution. However, the re- 
sults of [7f include ultrarelativistic near-threshold solu- 
tions for r = 1.99. This case does allow for the presence 
of the sonic point, so it may be that sonic-point kinks 
were not seen in this instance because p never became 
sufficiently steep at r s to excite the kink mode. 

In Fig. [TTJ we plot fits of the Type II scaling behavior 
of stars with three different values of T. The T = 2 
distribution is our fiducial system, while the other two — 
r = 1.88, 1.90 — use a different progenitor TOV solution 
having p c — 0.013 instead of p c — 0.05. The change 
in initial state was necessary since use of p c = 0.05 for 
r = 1.88 and 1.89 would not necessarily lead to Type II 
behavior using velocity induced collapse. In addition, the 
values 1.88 and 1.89 bound by a difference of ~ 1% the 
critical value, r c ~ 1.889, above which the kink mode 
becomes unstable [HJ . The scaling exponents derived 
from the two new sets of computations are 7(1.9) = 0.83 
and 7(1.88) = 0.81. The T = 1.9 scaling exponent is the 
same as that calculated with the ultrarelativistic PDEs 
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FIG. 17: Scaling behavior in T max for different values of F. 
The "Original" (solid line with circles) was made from runs 
with p c = 0.05 and T = 2 as before. Both the T = 1.88 
(dashed line with triangles) and V = 1.90 (dotted line with 
squares) runs used p c = 0.013 with {N a , Nt, N c , Ar a } = 
{420,700,20,0.005}. The scaling exponents are given in the 
text. 



in Q, while the T = 1.88 scaling exponent is consistent 
with the value obtained for T — 1.888, the value of T 
closest to 1.88 for which 7 was computed in Q- 

We find that there is no consistent behavior in code 
stability as the T c ~ 1.889 value is crossed. In fact, we 
find the opposite of the expected behavior: we are able 
to tune the kink-unstable (r = 1.90) data closer to p* 
than the kink-stable (r = 1.88) data. We therefore find 
it unlikely that the kink mode is the cause of our sonic 
point instability. 

Even if the kink mode is unstable for our EOS, tuning 
toward the CSS solution while in the presence of another 
unstable mode is not without precedent. For example, in 
a study of the spherically-symmetric general relativistic 
harmonic-map (nonlinear sigma model) , Liebling [53j | dis- 
covered that by judicious choice of an initial data family, 
he could tune to a critical solution that had been shown 
to have two unstable modes. However, it is unclear what 
relation this work might have to our current study, since 
the type of initial data that we are studying does not 
seem to have been chosen in any particularly special way 
(i.e. we suspect that the initial data that we have used 
is generic, whereas that used by Liebling to tune to the 
two- mode- unstable solution was, by construction, non- 
generic) . What is clear is that this issue requires further 
investigation, but we will leave that to future studies. 



V. CONCLUSION 

In this work, we simulated spherically-symmetric rel- 
ativistic perfect fluid flow in the strong-field regime of 
general relatively. Specifically, a perfect fluid that admits 
a length scale, for example one that follows a relativis- 
tic ideal-gas law, was used to investigate the dynamics 
of compact, stellar objects. These stars were modeled as 
neutron stars by using a stiff equation of state, approx- 
imating the behavior of some realistic state equations. 
Our models were then used to study the dynamics of neu- 
tron stars so far out of equilibrium that they are driven 
to gravitational collapse. 

Since these systems entail highly-relativistic fluid 
motions and strong, nonlinear effects from the fluid- 
gravitational interaction, a numerical treatment is chal- 
lenging. To achieve stable evolutions in near-luminal 
flows, while using high resolution shock capturing tech- 
niques, the primitive variable solver required improve- 
ments. In addition, an instability was found to develop 
in calculations near the threshold of black hole forma- 
tion, and this necessitated the use of new computational 
methods, that were only partially successful in stabilizing 
the calculations. 

We find our value for the scaling exponent, 7, given in 
Eq. ([57)1 agrees well with those found in [5J, and agrees 
with the value of 7 computed in [7| to within the uncer- 
tainty quoted in that work. We note that a discrepancy 
in the values of 7 computed using ideal-gas and ultrarela- 
tivistic fluids was observed in 0] , and this is also the case 
for our calculations. Our ultrarelativistic value, 7 = 0.97, 
agrees well with the value calculated in [7J], but deviates 
by an estimated 3 standard deviations from the value ex- 
tracted from our ideal-gas calculations. It is somewhat 
interesting, yet probably coincidental, that our results 
from the ideal-gas system of equations lead to estimates 
of 7 that agree with the perturbation calculations better 
than those values found from the ultrarelativistic PDE 
calculations. 

Our findings thus do not support some of the results 
found, and claims made by Novak 19] for the case of fluid 
collapse with an ideal-gas EOS and r = 2. This previ- 
ous work suggested that the Type II behavior observed 
in such a case was not well approximated by a univer- 
sal (with respect to initial data) ultrarelativistic limit. 
However, using different stars and velocity profiles, and 
by varying other aspects of the numerical model, we have 
found scaling behavior that is insensitive to approxima- 
tions made in the numerical solution, and which does 
appear to be universal with respect to families of initial 
data. Moreover, we have found that the scaling expo- 
nent and critical solution for the collapse governed by 
the ideal-gas EOS agrees well with their ultrarelativistic 
counterparts. 

Ultimately, it is our goal to expand the model a great 
deal, making the matter description more realistic and 
eliminating symmetry. As a first step, we wish to de- 
velop adaptive mesh refinement procedures for conserva- 
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tive systems that will be required to study critical phe- 
nomena of stellar objects in axial-symmetry [54j . 

It remains to be seen whether the universal scaling 
behavior we have observed is also seen with more real- 
istic state equations such as the one Novak used. Since 
accurate measurements of 7 have only been found for 
equations of state with constant adiabatic index T, and 
since 7 seems to only depend on T for perfect fluids, 
it will be interesting to investigate in detail what the 
scaling behavior — if any — will be like for realistic state 
equations with variable T. However, to the extent that 
any given realistic EOS admits a unique ultrarelativistic 
limit, characterized by a single value of T, we can expect 
to see universal Type II behavior of the sort discussed in 
this paper, for at least some collapse scenarios. 
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APPENDIX A: FINITE DIFFERENCE 
EQUATIONS AND APPROXIMATE RIEMANN 
SOLVERS 



We employ High-Resolution Shock Capturing (HRSC) 
algorithms to solve the equations of motion for the fluid 
([II]) . Such methods have become increasingly popular 
in the field of relativistic hydrodynamics since they are 
flux conservative, and ensure that discontinuities are well 
resolved and propagate at the correct speeds in the con- 
tinuum limit. A key ingredient to these schemes is their 
use of solvers for the Riemann problem at every cell in- 
terface. This is crucial for the conservative nature of 
these schemes since the solution to the Riemann prob- 
lem is always a weak solution of the hyperbolic conser- 
vation laws. The "high-resolution" aspect of the algo- 
rithms denotes that in regions where the grid functions 
are smooth, the integration procedure is at least 0(Ar 2 ) 
accurate. Many of the HRSC methods used in t his p aper 
have been used in previous works such as [26| . [l9( . and 
[25| to name only a few relevant sources. Also, excellent 
references on conservative methods for general systems 
of hyperbolic conservation equations have been written 
by LeVeque [H,|46|. 

Unlike finite difference methods, finite volume or con- 
servative methods calculate the cell-averages of grid func- 
tions instead of the grid functions themselves. The dif- 
ference equations for conservative methods are not de- 
rived from Taylor-series approximations of derivatives, 
but from differences of integrals. For instance, we differ- 



ence the fluid EOM (| 14[) in the following manner: 
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At first glance, this equation seems no different than a 
finite difference approximation of Eq. (I14|) . However, the 
difference between the two approaches becomes appar- 
ent when we examine how specific quantities in (|A1[) are 
defined, qf is the spatial average of q(r, t) over the cell 
centered at (r^, t"), is the spatio-temporal average of 
the source function tp centered at (r.i, i™ +1 / 2 ), and the 
numerical flux F™ +] y 2 is the time average of f(r 1+1 / 2 ,t) 
from f" to t n+1 . In practice, we approximate i/>™ as the 
source of the averages, ■0(q( r i J t™ +1//2 )). 

The techniques for determining the numerical flux are 
especially important since they set the spatial accuracy of 
the overall scheme and are primarily responsible for how 
well shocks are resolved by the method. Unless otherwise 
stated, all results presented in this paper were produced 
using an approximate Roe- type solver as outlined in 26 1 
for calculating the numerical flux. The Roe solver 35 1 
approximately solves the Riemann problem at each cell 
interface by casting the conservation equation (|14| into 
quasi-linear form. We approximate the Roe matrix A as 



with 



A(qW)=f; 



<H(q L + q") 



(A2) 



(A3) 



and where q L and q^ will be defined shortly. After solv- 
ing the linear Riemann problem, the numerical flux of is 
computed using the following expression (45| : 



F fc+ i/ 2 (i) = \ [f(q*+i/ 2 (t))+f(q£-i/2(*)) 



I ^m»7« 



(A4) 



Here, A m and rj m are the eigenvalues and right eigenvec- 
tors, respectively, of A, q L and q fl are, respectively, the 
values of q to the left and right of the cell boundary, and 
Lj m are the decomposed values of the jumps in the space 
of characteristic values: 



- q L = Y2 u ™rin 



(A5) 



In order to calculate all quantities associated with A, 
such as A m and r] m , we use the average of the left and 
right states, q, defined by (|A3|) . 
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As far as we know, the eigenvectors for our new for- 
mulation, ([14]) and (jTTJ) , have not previously been pub- 
lished, so for completeness we present them here. Since 
the transformation from {D, S, t} to {D, II, $} is linear, 
the eigenvalues remain the same: 



Ai = v 



v±c s 

\2 = A± = — — 

3 1 ± V C, 



(A6) 



Using Maple and assuming the the ideal-gas EOS ([3]), 
we have calculated the left and right eigenvectors; the 
Marquina flux formula requires the right eigenvectors 
36] . Using the typical normalization for the eigenvectors 

(2) 

(rim = A m ), leads to a very complicated set of eigenvec- 
tors. Hence, we used the normalizations 



V 



(i) _ 



1 



V m . 



(A7) 



which leads to significant simplification. The right eigen- 
vectors become: 



r\i = r] ± 



1 

W(l+v) _ - 
a 

W(l-v) _ - 
a 

1 

2^Ml=Fc.)-l 



(A8) 



(A9) 



while the associated left eigenvectors are 

1 



li 



2h c 



2h t 



kW(l-v) 



kW(l + v) 



,(A10) 



1± 



2// c- 



aW (k =p vc s ) — k 
\aW{l-v) (k±c s ) 
\aW{l + v) (kTc s ) 
In the above expressions we have 

(r - 1)TP 



(All) 



ct = 



he 



(r — i)p + rp 
r-i, 

TP 

Po 



(A12) 



The simple form the numerical flux (|A4[) takes in the 
approximate Roe method allows us to subtly alter the 
equations of motion in a way that greatly improves the 
regularity of the conserved variables near the origin [25j . 
Let us first note that the flux term in Eq. (I14|) can be 
expanded in a manner that yields a new EOM: 



if, (A13) 



where q remains is unchanged from Eq. (|17[) . and 
f(i) = 



f(2) 



4> 



Dv 

v{n + p) 

v($ + P) 


p 

-p 



e 

-e 



(A14) 



This new formulation eliminates the inexact cancellation 
of the 2PX/r terms from the flux and the source that 
would normally arise from truncation error. Note that 
the eigensystem, used by Roe's numerical flux function 
(|A4|) . is still calculated from the total flux function, f = 

f(l) +f(2). 

Due to the finite precision of the calculations and the 
nature of the numerical methods employed, the evacua- 
tion of fluid often "overshoots" the vacuum state gener- 
ating negative pressures or densities, which in turn leads 
to problems such as a complex c s or super-luminal ve- 
locities. In order to alleviate such problems we require 
the dynamic fluid quantities — the conserved variables q 
(fT7|) — to have values greater than or equal to a so-called 
"floor" state. In order to determine the floor state, we 
require P. p Q > and \v\ < 1 which implies that 



D , (r±\S\) > 



(A15) 



Using the transformed ("new") variables II, $, we imple- 
ment this requirement in the following way 

D = max (D,S) , (A16) 
II = max (n + D, 25) - D , (A17) 
$ = max ($ + D, 25) - D , (A18) 

where 5 is the adjustable floor parameter. Notice that II 
and $ need not remain positive since r < is physical 
as long as E > 0. 



APPENDIX B: BOUNDARY CONDITIONS 

For the outer boundary condition of our fluid quanti- 
ties, we use the typical outflow condition where the fluid 
quantities associated with the last physical cell are copied 
into so-called ghost cells (i.e. first order extrapolation). 
Our experience, as well as that of others, indicates that 
this condition is fairly robust and non-reflective so no 
other methods were tested or used. 

The regularity conditions at the origin are, however, 
more sophisticated. Since our fluid grid functions are 
defined with respect to a grid that is offset from the ori- 
gin, typical 0(Ar 2 ) regularity conditions are not as well- 
behaved as they are for origin-centered cells. Hence, we 
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have found it helpful to use higher-order, conservative in- 
terpolation for the fields on the first physical cell. Since 
the fluid fields, q^, are to be interpreted as cell-averages 
of conserved functions, which we will call Q(j~), an in- 
terpolation is said to be conservative if the integral of a 
function on a local domain is conserved by the interpo- 
lation procedure. We first assume that the interpolation 
function Qj(r) that is associated with a cell Ci has a 
polynomial expansion of degree N — 1: 



N-l 



(Bl) 



with N coefficients a„. These coefficients are found by 
demanding that Q; maintains conservation locally. That 
is, a set Si of TV cells is chosen in the neighborhood of cell 
Ci, and we require that reproduces the known values 
q/c, where Ck £ Si. Specifically, the coefficients a„ are 
calculated by solving the following set of N equations: 



q_k = ^ [ Qi(r) dV 
Vk J Vk 



iy»3 /v>3 

'fe+l/2 1 k-1/2 



(B2) 



N-l 

x ^ a n 

n=0 



r k + l/2 



(r — 7*j) r dr 



Th-X/2 



for all Ck € Si . Since this interpolation procedure is used 
at the origin where local flatness is demanded, we can 
assume a(r, t) = 1 for r as with negligible effect (in 
principle, a(r, t) should appear in the above integral as 
part of the volume element). Once Eq. (|B2j) is solved for 
the coefficients, a n , the interpolation procedure is com- 



pleted by using ()B2|) to determine values q 3 for a cell 

Cj i Si. 

From the demand of regularity at the origin, the fields 
p 0l P,D,r are all even in r as r — ► 0, while v and S 
are odd. Thus, a n = for odd n in the interpolation 
function of the even fields, and a n = for even n in the 
odd interpolations. In our case, the cells lying nearest 
the origin are spaced uniformly and we use N — 4. For 
even functions the boundary condition then becomes 



3311 q 2 - 2413 q 3 + 851 q 4 - 122 q 5 



1627 

while for odd functions we have 

35819 q 2 - 16777 q 3 + 4329 q 4 



q 5 



36883 



(B3) 



(B4) 



Here, cell C\ is the innermost cell, located at r% = Ar/2. 
Since II and $ are combinations of even and odd func- 
tions, their regularity conditions are not as straightfor- 
ward to compute. To determine their behavior at the 
origin, we first calculate the interpolated values of r and 
S at C\, since their regularity behavior is known. Then, 
n and $ are calculated for C\ from the definitions CE 
using the interpolated values of r and S. 



In contrast to the fluid variables, the metric grid func- 
tions are defined on a grid centred on the origin. This 
enables us to use straightforward discrete expressions to 
ensure the regularity of a and a. In solving the Hamilto- 
nian equation (|20[) , we demand that spacetime be locally 
flat at the origin; this implies a(Q,t) = 1. The slicing 
condition (|2ip is solved by integrating inward from the 
outer boundary, making use of the freedom we have in 
relabeling constant t surfaces. If we assume that all the 
matter remains within our grid, then the metric exterior 
to the grid is a piece of the Schwarzschild solution. Since 
the Schwarzschild metric is asymptotically flat, we can 
rescale a after the solution of the slicing equation so that 
our metric is equivalent to the standard Schwarzschild 
form at r max . Thus, our boundary condition for a is 



C^(^max) 



1 



a(r n 



(B5) 



The physical meaning of this commonly-adopted param- 
eterization is that coordinate time and proper time coin- 
cide as r — > oo. 



APPENDIX C: NUMERICAL TESTS 

Here, we present a series of tests that verify that our 
code solves the equations we claim that it solves, and 
that our discrete solutions converge as expected in the 
continuum limit. 




FIG. 18: Riemann solution using the approximate Roe 
method with initial data {P, v,p }(x < 0.5) = {100,0,1}, 
{P,v, Po }(x > 0.5) = {1,0,1}, T = 5/3, using 200 cells. 
P(x)/120 (circles), p a (x)/6 (triangles), and v(x) (x) are plot- 
ted at t = 0.4. The lines correspond to the exact solution. 

In Fig. [TS1 we show solutions of a Riemann problem 
generated using 1) the approximate Roe solver described 
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previously, and 2) Marquina's method. Here, the Rie- 
mann problem was set up in the middle of the grid, 
i.e. at x = 0.5. The solid line shows the exact solution 
of the Riemann problem calculated by a routine given 
in [49(. The approximate solutions compare favorably 
to the exact solution, especially in smooth regions where 
the discrete solutions should be close to the exact solu- 
tion. It seems that Roe's method produces a slight Gibbs 
phenomenon near the origin of the rarefaction fan, while 
Marquina's method results in undershooting in that re- 
gion. However, the two methods produce nearly identical 
results near the shock and contact discontinuity. 

To illustrate convergence properties of our discrete ap- 
proximations, we show the results of a convergence test 
of our code in Figs. |2"0H2"21 for the quantities D, n, and 
$ (fl7|) , respectively. The scaled error estimates shown 
in the top panels of each figure demonstrate how the 
computed value of each dynamical variable exhibits the 
expected dependence on the fundamental grid spacing, 
h. Specifically, so long as the flow is smooth (including 
the initial conditions) , and our spatial and temporal grid 
spacings are characterized locally by a single discretiza- 
tion scale, h, then because our scheme is second order, 
we expect Richardson expansions of the form: 

lim f h (t, r) = f(t, r) + h 2 e 2 (t, r) + ■ ■ ■ (CI) 

where f(t,r) represents any dynamical fluid variable, 
f h (t,r) is the discrete approximation to that variable, 
and e2(t,r) is an /i-independent function with smooth- 
ness comparable to /. The data shown in the plot 
have been extracted at a time before any discontinuities 



FIG. 20: Convergence test for the fluid variable D. The top 
panel shows In | D 8h - D 4h \ (solid line), In |4 (D 4h — D 2h ) j 
(dots), In j 16 (D 2h - D h ) | (dashes) which have been scaled so 
that they will coincide if they are well-described by Richard- 
son expansions of the form (|C1|I . Quantities with super- 
script "Ih" have been calculated using grid point separations 
I times larger than the fiducial I — 1 quantities. The bot- 
tom panel shows D(r, 0) (dashes) and D(r,t) (solid line), 
where t is the time at which we performed the convergence 
test. The initial data consisted of a self-gravitating fluid 
specified by a Gaussian function for p a centered at r — 0.1 
with an initial linear velocity profile. The initial grid used 
for the coarsest solution shown is defined by the parameters 
{N a , N b , N c , Ar a } = {200, 300, 20, 0.005}. 



regridding procedure described above, we performed the 
convergence test at a time following the first grid refine- 
ment. From these results, it is evident that our numerical 
methods are 2nd-order accurate for smooth flows. 



APPENDIX D: PSEUDO-CODE FOR THE 
PRIMITIVE VARIABLE CALCULATION 

The primitive variables are calculated from the 
conserved variables using a one-dimensional Newton- 
Raphson method that locates the value of H that mini- 
mizes the residual given in Eq. (|28p . The equations for 
calculating the primitive variables and other quantities 
as a function of the value of H in a given iterations are 
listed in Table HVl for different limits of A. 
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FIG. 21: Convergence test for the fluid variable II. The top 
panel shows the scaled error estimates described in the cap- 
tion of Fig. l20l The bottom panel shows II(r, 0) (dashed) and 
II(r, t) (solid), where t is the time at which convergence is 
tested. 



0.0005 




FIG. 22: Convergence test for the fluid variable The top 
panel shows the scaled error estimates described in the cap- 
tion of Fig. [20] The bottom panel shows $(r, 0) (dashed) and 
$(r, t) (solid), where t is the time at which convergence is 
tested. 
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